#include "mfAddFuncCalcTask.h"
#include "kernel/feMaterial.h"
#include "kernel/fe1dParabolicTask.h"
#include "kernel/feBordCond.h"
#include "kernel/fePolylineGenerator.h"
#include "ConductiveLayer.h"
#include "kernel/fePolyline.h"
#include "kernel/feObserver.h"
#include "kernel/feMeshFunction.h"

#include <vector>
#include <math.h>


mfAddFuncCalcTask::mfAddFuncCalcTask()
{
	Layers = NULL;
	DistMesh = NULL;
	DistFunc = NULL;
	NumLayers = -1;
	Sigma0 = -1;
	MaxZ = 0;
}


//------------------------------------------------------------------------------------------------------------------------------------------
void mfAddFuncCalcTask::CalcSigma0()
{
	//check num layers
	if(NumLayers == 1) 
	{
		Sigma0 = Layers[0].Sigma;
		MaxZ = 3000/M;
		return;
	}


	double firsthe = 3.0/M;
	double acche = 0;
	

	//calc first sigma as average through the layers
	Sigma0 = 0;
	for(int i=0; i<NumLayers; i++)
	{
		Sigma0 += Layers[i].Sigma;
	}
	Sigma0 /= NumLayers;


	//scale layers and calc Ne, He for Zero Sigma
	int ne = -1;
	double dh = 0;
	for(int i=0; i<NumLayers; i++)
	{
		if(Layers[i].H < MF_INFINITY_H)
		{
			firsthe -= Layers[i].H*Layers[i].Sigma/Sigma0;
			if(firsthe < 0)
			{
				firsthe += Layers[i].H*Layers[i].Sigma/Sigma0;
				dh = firsthe*Sigma0/Layers[i].Sigma;
				acche += dh;
				ne = i;
				break;
			}

			acche += Layers[i].H;
		}
		else
		{
			dh = firsthe*Sigma0/Layers[i].Sigma;
			acche += dh;
			ne = i;
			break;
		}		
	}	
	
	if(ne == -1) 
		ne = NumLayers;

	
	//save maxz
	MaxZ = acche*100;



	double maxsigma = 0;
	double minsigma = 1e48;
	for(int i = 0; i <= ne; i++)
	{
		if(maxsigma < Layers[i].Sigma)
			maxsigma = Layers[i].Sigma;
		if(minsigma > Layers[i].Sigma)
			minsigma = Layers[i].Sigma;
	}


	//try to eliminate big contrast
	if(MaxZ < 10*Layers[0].H && maxsigma / minsigma > 10 )
	{
		for(int i = 0; i <= ne; i++)
		{
			if(maxsigma > 10*Layers[i].Sigma)
			{
				Sigma0 = Layers[i].Sigma;				
				return;
			}
		}
	}
	

	//calc zero sigma as layers' power average
	if(dh != 0)
		Sigma0 = dh*Layers[ne].Sigma;
	for(int i = 0; i < ne; i++)
		Sigma0 += Layers[i].Sigma*Layers[i].H;
	Sigma0 /= acche;	
}


//------------------------------------------------------------------------------------------------------------------------------------------
bool mfAddFuncCalcTask::Load()
{
	if(!Helper->Load()) return false;
	return true;
}


//------------------------------------------------------------------------------------------------------------------------------------------
bool mfAddFuncCalcTask::Init(ConductiveLayer *iLayers, int iLayersCount, double iTimeCoeff, double iFirstTimeStep, 
					   int iLayersTime, double iM, TaskHelper *iHelper, bool iUseMainField)
{
	//copy coeffs
	Helper = iHelper;
	M = iM;
	FirstTimeStep = iFirstTimeStep;
	TimeCoeff = iTimeCoeff;
	LayersTime = iLayersTime;
	UseMainField = iUseMainField;

	NumLayers = iLayersCount;
	Layers = iLayers;
	for(int i=0; i<NumLayers; i++)
	{
		if(Layers[i].H >= MF_INFINITY_H)
			Layers[i].H = 1e20;
	}
	

	//calc zero sigma
	CalcSigma0();


 	//check maxz
	if(UseMainField && MaxZ < Layers[0].H)
	{
		NumLayers = 1;
		Sigma0 = Layers[0].Sigma;
	}

	
	
	//load main field
	if(UseMainField)
	{
		if(!Load())
			UseMainField = false;
	}


	//create mesh
	CreateMesh();
	Helper->SetMesh(DistMesh);
	

	//done
	return true;
}


//------------------------------------------------------------------------------------------------------------------------------------------
void mfAddFuncCalcTask::CreateMesh()
{
	if(UseMainField && NumLayers == 1) return;

// 	if(UseMainField)
// 	{
// 		GenAdoptedMesh();
// 		return;
// 	}
	
	double stepcoef;
	double fstep;
	if(UseMainField)
	{
		stepcoef = 2.0;
		fstep = 0.1 > Layers[0].H ? Layers[0].H : 0.1;
	//	double npoints = 20;
	//	fstep = MaxZ*(1-stepcoef)/(1-pow(stepcoef, npoints));
	}
	else
	{
		stepcoef = 1.15;
		fstep = 0.1;
	}


	//check mesh
	if(!DistMesh)
		DistMesh = new mfMesh();


	//clear task
	DistMesh->ClearMesh();



	//generate mesh
	int pcount = 1;
	int elcount = 0;
	double z = 0;
	DistMesh->PushNode(&fePoint(z));
	double step = fstep;
	double scale = 1;
	for(int i=0; fabs(z-MaxZ)>1e-5*MaxZ; i++)
	{
		//create and add material
		feSharedPtr<feMaterial> m;
		scale = Helper->GetScaleSigma(i);
		if(UseMainField)
		{

			m = new feMaterial(scale, M*M*scale, MU0*Layers[i].Sigma*scale, Helper->GetRightFunction(), Helper->GetScaleFCoeff(i, Sigma0), M, M*M/Sigma0);
		}
		else
		{
			m = new feMaterial(scale, M*M*scale, MU0 * Layers[i].Sigma*scale, 0, 0, 0);
		}

		double h = min(Layers[i].H, MaxZ-z);		


		//gen mesh for the layer
		double z0 = z;	
		while(fabs(z-z0-h) > 1e-8 && fabs(z-MaxZ)>1e-5*MaxZ)
		{
			double st = step;
			if(z - z0 + step > h)
				st = h - z + z0;

			//add nodes
			DistMesh->PushNode(&fePoint(z+st/2));
			
			z += step;

			//check layer's border
			if(z-z0 > h)
				z = z0 + h;				

			DistMesh->PushNode(&fePoint(z));

			//add element
			feQuadraticElement1d e(pcount-1, DistMesh, m);
			e.SetRightCoeff(Helper->GetStiffElemCoeff(), Helper->GetWeightElemCoeff());
			DistMesh->PushElement(&e);

			
			//update vars			
			elcount++;
			pcount+=2;

			step = fstep*pow(stepcoef, elcount);
		}

// 		if(z != MaxZ && Helper->NeedDoubleNodes())
// 		{
// 			//double nodes on layer border
// 			//DistMesh->PushNode(&fePoint(z));
// 			jump_nodes[jnnum] = pcount-1;
// 			jnnum++;
// 			//pcount++;
// 		}
	}
}


//------------------------------------------------------------------------------------------------------------------------------------------
bool mfAddFuncCalcTask::Solve()
{
	if(!DistFunc)
		DistFunc = new feSimple1dMeshFunction(DistMesh);
	else
		DistFunc->Clear();

	//don't need to solve anything
	if(UseMainField && NumLayers == 1) return true;


	//get X0
	X0 = Helper->GetX0(UseMainField);

	
	//create border condition
	feBordCond **bcond;	
	int bcnum;
	bcond = Helper->GetBcs(UseMainField, bcnum);
	


	fe1dParabolicTask task;
	//init sover
	if(!UseMainField)
	{
		if(!task.Init(DistMesh, bcond, bcnum, X0, FirstTimeStep, LayersTime, TimeCoeff))
			return false;
	}
	else
	{
		//get layers count and layers pass
		double tc = M*M/Sigma0;
		int lt = 1;
		double maxt = FirstTimeStep * (1-pow(TimeCoeff, LayersTime))/(1-TimeCoeff);
		for(int i=0; i<Helper->GetRightFunction()->LayersCount(); i++)
		{
			if(Helper->GetRightFunction()->GetTime(i) > maxt*tc)
			{
				lt = i/LayersTime; //assuming that we will use every lt's layer in time mesh
				break;
			}
		}
		if(lt == 0) lt = 1;
						
		
		if(!task.Init(DistMesh, bcond, 2, X0, Helper->GetRightFunction(), lt, LayersTime, M*M/Sigma0))
			return false;
	}


	//solve
	if(!task.SolveNonStationary(DistFunc))
		return false;

	return true;
}


//------------------------------------------------------------------------------------------------------------------------------------------
double mfAddFuncCalcTask::GetDuDt(double iZ, double iTime)
{ 	
 	if(UseMainField && NumLayers == 1)
 		return M*M/Sigma0*Helper->GetMfDuDtFunction()->GetU(&fePoint(iZ*M), iTime*M*M/Sigma0);
  	else if(UseMainField)
		return DistFunc->GetDuDt(&fePoint(iZ), iTime) + M*M/Sigma0*Helper->GetMfDuDtFunction()->GetU(&fePoint(iZ*M), iTime*M*M/Sigma0);
	else 
		return DistFunc->GetDuDt(&fePoint(iZ), iTime);
}


//------------------------------------------------------------------------------------------------------------------------------------------
double mfAddFuncCalcTask::GetU(double iZ, double iTime)
{
	if(UseMainField && NumLayers == 1)
		return Helper->GetMfUFunction()->GetU(&fePoint(iZ*M), iTime*M*M/Sigma0);
	else if(UseMainField)
		return DistFunc->GetU(&fePoint(iZ), iTime) + Helper->GetMfUFunction()->GetU(&fePoint(iZ*M), iTime*M*M/Sigma0);
	else 
		return DistFunc->GetU(&fePoint(iZ), iTime);
}


//------------------------------------------------------------------------------------------------------------------------------------------
double mfAddFuncCalcTask::GetDuDz(double iZ, double iTime)
{	
	if(UseMainField && NumLayers == 1)
		return M*Helper->GetMfUFunction()->GetDuDx(&fePoint(iZ*M), iTime*M*M/Sigma0);
	else if(UseMainField)
		return DistFunc->GetDuDx(&fePoint(iZ), iTime) + M*Helper->GetMfUFunction()->GetDuDx(&fePoint(iZ*M), iTime*M*M/Sigma0);
	else 
		return DistFunc->GetDuDx(&fePoint(iZ), iTime);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
bool mfAddFuncCalcTask::CalcMainField()
{
	if(!ZTaskHelper::CalcMainField()) return false;
	if(!XTaskHelper::CalcMainField()) return false;
	return true;
}


//------------------------------------------------------------------------------------------------------------------------------------------
mfAddFuncCalcTask::~mfAddFuncCalcTask()
{
	if(DistFunc) delete DistFunc;
}